Note: This tutorial assumes that you have successfully installed ABSEIR, and are at least passingly familiar with compartmental models. Some introductory information is available in this vignette.
Data for this example were obtained from the measlesWeserEms data set from the surveillance package.
# load the ABSEIR library
library(ABSEIR)
## Loading required package: Rcpp
## Loading required package: parallel
## Loading required package: compiler
library(surveillance)
## Loading required package: sp
## Loading required package: xtable
## Loading required package: polyCub
## This is surveillance 1.13.0. For overview type 'help(surveillance)'.
# Read in data
data(measlesWeserEms)
# Identify cases
cases<-measlesWeserEms@observed
epidemic.start = min(which(apply(cases, 1, max) > 0))
cases = cases[(epidemic.start-1):nrow(cases),]
# Obtain distance matrix
neighbourhood<-measlesWeserEms@neighbourhood
# Week
week <- measlesWeserEms@epoch[(epidemic.start-1):length(measlesWeserEms@epoch)]
# Vaccination and population data set
vaccine.data <- measlesWeserEms@map@data
vaccine.data$adminID <- rownames(vaccine.data)
# Population size
N <- vaccine.data$POPULATION
# Check that spatial unit ordering makes sense
if (!all(vaccine.data$adminID == colnames(neighbourhood))){
stop("Error, make sure spatial unit ordering is consistent.")
}
# Make cumulative
weser.data_model = DataModel(Y=apply(cases, 2,cumsum),
type = "identity", # Assume data is correct
compartment = "I_star", # Data related to new infections
cumulative = TRUE # Not reported on cumulative scale
)
I0 = (apply(cases[1:3,], 2, max) > 0)*2
E0 = I0
R0 = 0*I0
S0 = N-E0-I0-R0
weser.initial_values = InitialValueContainer(S0 = S0,
E0 = E0,
I0 = I0,
R0 = R0)
# No reinfection
weser.reinfection_model = ReinfectionModel("SEIR")
To specify the exposure process, we need to create the exposure design matrix. This matrix of temporal and spatial covariates must have a particular structure. Each of the \(J\) spatial locations receives a \(T\times p\) design matrix, where \(p\) is the number of covariates. These matrices are then ‘stacked’ row-wise in the same spatial ordering as presented in the columns and rows of the neighborhood matrix (defined above).
Here, we include the following covariates:
For more information on these covariates, see the documentation for the ‘measlesWeserEMS’ data.
As you can see, while ABSEIR doesn’t directly fit SEIRV models (which incorporate a vaccination compartment), we can include vaccination data into the exposure probability structure, effectively giving vaccinated persons a different probability of contracting the pathogen of interest.
n.locations <- ncol(neighbourhood)
n.timepoints <- length(week)
time_invariant_covariates <- data.frame(intercept = rep(1, nrow(vaccine.data)),
popDensity = vaccine.data$POPULATION/vaccine.data$AREA,
proportionVaccineCard = vaccine.data$vaccdoc.2004,
proportionVaccine1 = vaccine.data$vacc1.2004,
proportionVaccine2 = vaccine.data$vacc2.2004)
time_varying_covariates <- data.frame(sin_component = sin(week/52*2*pi),
cos_component = cos(week/52*2*pi),
trig_interact = sin(week/52*2*pi)*cos(week/52*2*pi))
exposure.design.matrix <- as.matrix(
cbind(
time_invariant_covariates[rep(1:n.locations, each = n.timepoints),],
time_varying_covariates[rep(1:n.timepoints, n.locations),]
)
)
## Build the exposure model
weser.exposure_model <- ExposureModel(X = exposure.design.matrix,
nTpt = n.timepoints,
nLoc = n.locations,
betaPriorPrecision = rep(1, ncol(exposure.design.matrix)),
betaPriorMean = rep(0, ncol(exposure.design.matrix)))
## [1] "Assuming equally spaced count data."
We now need to specify the contact structures. Here, we’ll build both a gravity model (which depends on population size and distance), and a simple CAR model (which just depends on shared borders between regions).
# Build a gravity model, in which the contact process between a pair of
# spatial locations is proportional to the product of their populations divided
# by the squared 'distance'
pop.matrix <- matrix(N, nrow = length(N), ncol = length(N))
gravityModel <- (pop.matrix * t(pop.matrix))/neighbourhood^2
diag(gravityModel) <- 1
# Rescale
maxRowSum <- max(apply(gravityModel,1,sum))
gravityModel <- gravityModel/maxRowSum
weser.distance_model <- DistanceModel(list(gravityModel),
priorAlpha = 1,
priorBeta = 1)
# Build a simpler contact model, in which the contact probability between
# a pair of spatial locations is only nonzero when the distance is equal to 1
# (CAR specification)
weser.CAR_model <- DistanceModel(list((neighbourhood == 1)*1),
priorAlpha = 1,
priorBeta = 1
)
# 9-12 day latent period
# Infectious
weser.transition_priors = ExponentialTransitionPriors(p_ei = 0.8, # Guess at E to I transition probability (per week)
p_ir = 0.8, # Guess at I to R transition probability (per week)
p_ei_ess = 10, # confidence
p_ir_ess = 10 # confidence
)
weser.sampling_control <- SamplingControl(seed=123124,
n_cores = 14,
algorithm = "Beaumont2009",
params = list(
batch_size = 2000,
init_batch_size = 1000000,
epochs = 1e6,
shrinkage = 0.99,
max_batches = 200,
multivariate_perturbation=FALSE
))
# Cache the model objects so we don't need to re-run the analysis for things like format tweaks
if (!file.exists( "weserModel1.rda")){
t1 <- system.time(
{
weser.model1 <- SpatialSEIRModel(data_model = weser.data_model,
exposure_model = weser.exposure_model,
reinfection_model = weser.reinfection_model,
distance_model = weser.distance_model,
transition_priors = weser.transition_priors,
initial_value_container = weser.initial_values,
sampling_control = weser.sampling_control,
samples = 50,
verbose = FALSE)
}
)
save("weser.model1", "t1", file = "weserModel1.rda")
} else {
load("weserModel1.rda")
}
if (!file.exists( "weserModel2.rda")){
t2 <- system.time(
{
weser.model2 <- SpatialSEIRModel(data_model = weser.data_model,
exposure_model = weser.exposure_model,
reinfection_model = weser.reinfection_model,
distance_model = weser.CAR_model,
transition_priors = weser.transition_priors,
initial_value_container = weser.initial_values,
sampling_control = weser.sampling_control,
samples = 50,
verbose = 2)
}
)
save("weser.model2", "t2", file = "weserModel2.rda")
} else {
load("weserModel2.rda")
}
print(t1)
## user system elapsed
## 10504.509 16.417 760.814
print(t2)
## user system elapsed
## 23474.217 34.397 1701.912
summary(weser.model1)
## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Exposure Process Parameters: 8
## Reinfection Model Parameters: 0
## Spatial Parameters: 1
## Transition Parameters: 2
##
##
## Parameter Estimates:
## Mean SD 95% LB 95% UB
## Beta_SE_1 -0.350 0.719 -1.752 0.946
## Beta_SE_2 0.475 0.790 -1.299 1.995
## Beta_SE_3 -0.632 0.615 -1.944 0.497
## Beta_SE_4 0.108 0.518 -0.966 0.923
## Beta_SE_5 0.708 0.890 -1.123 2.057
## Beta_SE_6 0.050 0.203 -0.359 0.414
## Beta_SE_7 0.224 0.318 -0.466 0.830
## Beta_SE_8 0.124 1.570 -3.647 2.266
## rho_1 0.116 0.036 0.043 0.183
## gamma_EI 1.275 0.660 0.244 2.341
## gamma_IR 1.574 0.534 0.757 2.575
summary(weser.model2)
## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Exposure Process Parameters: 8
## Reinfection Model Parameters: 0
## Spatial Parameters: 1
## Transition Parameters: 2
##
##
## Parameter Estimates:
## Mean SD 95% LB 95% UB
## Beta_SE_1 -0.201 0.307 -0.792 0.495
## Beta_SE_2 0.103 1.181 -2.181 2.461
## Beta_SE_3 -0.176 0.386 -0.980 0.509
## Beta_SE_4 -1.018 0.578 -1.833 0.084
## Beta_SE_5 0.921 0.331 0.300 1.708
## Beta_SE_6 0.162 0.425 -0.648 0.883
## Beta_SE_7 -0.378 0.406 -1.008 0.349
## Beta_SE_8 -0.139 0.868 -1.772 1.392
## rho_1 0.026 0.010 0.009 0.044
## gamma_EI 1.595 0.396 0.728 2.179
## gamma_IR 1.280 0.754 0.417 2.709
This function will plot cases each of the spatial locations for a model, add a black line for the mean posterior predictive number of cases, and add blue dashed lines for the 99% credible interval for the number of cases.
makePlots <- function(modelObj, nm){
sims <- epidemic.simulations(modelObj, replicates = 50)
Is <- lapply(sims$simulationResults, function(x){x$I_star})
Is <- array(Reduce(c, Is), dim = c(nrow(Is[[1]]),
ncol(Is[[2]]),
length(Is)))
Ism <- apply(Is, 1:2, mean)
Islb <- apply(Is, 1:2, quantile, probs = c(0.025))
Isub <- apply(Is, 1:2, quantile, probs = c(0.975))
plotLocation <- function(x, model){
plot(cases[,x], ylim = c(0, max(Isub[,x])),
main = paste(model, ": Observed and Posterior Predictive Simulation.\n location ",
colnames(neighbourhood)[x], sep = ""))
lines(Ism[,x], col = rgb(0,0,0,0.8), lwd = 2)
lines(Islb[,x], col = rgb(0,0,0.5,0.8), lwd = 1, lty = 2)
lines(Isub[,x], col = rgb(0,0,0.5,0.8), lwd = 1, lty = 2)
#apply(Is, 3, function(i){
# lines(i[,x], col = rgb(0,0,0,0.1))
#})
points(cases[,x], pch = 16, col = "blue")
}
for (i in 1:ncol(neighbourhood)){
plotLocation(i, nm)
}
}
makePlots(weser.model1, "Dist")
makePlots(weser.model2, "CAR")
comps <- compareModels(modelList = list(weser.model1, weser.model2), n_samples = 100)
## [1] "Assuming equal prior probabilities."
rownames(comps) <- colnames(comps) <- c("Distance", "CAR")
print(comps)
## Distance CAR
## Distance 1.0 0.1315789
## CAR 7.6 1.0000000
It looks like the preferred model is the CAR model.